module mod_migrate
  ! finite difference step

  implicit none
  real(8)::dt,gamma

contains

  subroutine migrate()
    ! ch-->h-->ce-->e-->ch
    use mod_para
    use mod_workspace
    use mod_grid
    use mod_material
    implicit none
    integer::i,j,k
    do k=nz1,nz2
       do j=ny1,ny2
          do i=nx1,nx2
             hx(i,j,k)=g11(i,j,k)*chx(i,j,k)+ &
                  0.25D0*g12(i,j,k)*(chy(i-1,j+1,k)+chy(i,j+1,k)+chy(i-1,j,k)+chy(i,j,k))+ &
                  0.25D0*g13(i,j,k)*(chz(i-1,j,k+1)+chz(i,j,k+1)+chz(i-1,j,k)+chz(i,j,k))
             hy(i,j,k)=g22(i,j,k)*chy(i,j,k)+ &
                  0.25D0*g23(i,j,k)*(chz(i,j-1,k+1)+chz(i,j-1,k)+chz(i,j,k+1)+chz(i,j,k))+ &
                  0.25D0*g12(i,j,k)*(chx(i+1,j-1,k)+chx(i+1,j,k)+chx(i,j-1,k)+chx(i,j,k))
             hz(i,j,k)=g33(i,j,k)*chz(i,j,k)+ &
                  0.25D0*g13(i,j,k)*(chx(i+1,j,k-1)+chx(i+1,j,k)+chx(i,j,k-1)+chx(i,j,k))+ &
                  0.25D0*g23(i,j,k)*(chy(i,j+1,k-1)+chy(i,j+1,k)+chy(i,j,k-1)+chy(i,j,k))
          enddo
       enddo
    enddo

    do k=nz1,nz2
       do j=ny1,ny2
          do i=nx1,nx2
             cex(i,j,k)=2D0*dt*(hy(i,j,k-1)-hy(i,j,k)-hz(i,j-1,k)+hz(i,j,k)) / ( jac(i,j,k)*( 2*gamma+dt*sigma(i,j,k) ) ) &
                  +cex(i,j,k)*( 2D0*gamma-dt*sigma(i,j,k) )/( 2D0*gamma+dt*sigma(i,j,k) )
             cey(i,j,k)=2D0*dt*(-hx(i,j,-1+k)+hx(i,j,k)+hz(-1+i,j,k)-hz(i,j,k))/(jac(i,j,k)*(2*gamma + dt*sigma(i,j,k) ) ) &
                  +cey(i,j,k)*( 2D0*gamma-dt*sigma(i,j,k) )/( 2D0*gamma+dt*sigma(i,j,k) )
             cez(i,j,k)=2D0*dt*(hx(i,j-1,k)-hx(i,j,k)-hy(-1+i,j,k)+hy(i,j,k))/(jac(i,j,k)*(2*gamma + dt*sigma(i,j,k) ) ) &
                  +cez(i,j,k)*( 2D0*gamma-dt*sigma(i,j,k) )/( 2D0*gamma+dt*sigma(i,j,k) )
          enddo
       enddo
    enddo

    do k=nz1,nz2
       do j=ny1,ny2
          do i=nx1,nx2
             ex(i,j,k)=g11(i,j,k)*cex(i,j,k)+ &
                  0.25D0*g12(i,j,k)*(cey(i+1,j-1,k)+cey(i,j-1,k)+cey(i+1,j,k)+cey(i,j,k))+ &
                  0.25D0*g13(i,j,k)*(cez(i+1,j,k-1)+cez(i,j,k-1)+cez(i+1,j,k)+cez(i,j,k))
             ey(i,j,k)=g22(i,j,k)*cey(i,j,k)+ &
                  0.25D0*g23(i,j,k)*(cez(i,j+1,k-1)+cez(i,j+1,k)+cez(i,j,k-1)+cez(i,j,k))+ &
                  0.25D0*g12(i,j,k)*(cex(i-1,j+1,k)+cex(i-1,j,k)+cex(i,j+1,k)+cex(i,j,k))
             ez(i,j,k)=g33(i,j,k)*cez(i,j,k)+ &
                  0.25D0*g13(i,j,k)*(cex(i-1,j,k+1)+cex(i-1,j,k)+cex(i,j,k+1)+cex(i,j,k))+ &
                  0.25D0*g23(i,j,k)*(cey(i,j-1,k+1)+cey(i,j-1,k)+cey(i,j,k+1)+cey(i,j,k))
          enddo
       enddo
    enddo

    do k=nz1,nz2
       do j=ny1,ny2
          do i=nx1,nx2
             chx(i,j,k)=(ey(i,j,k+1)-ey(i,j,k)+ez(i,j,k)-ez(i,j+1,k))*dt/(mu*jac(i,j,k))+chx(i,j,k)
             chy(i,j,k)=(ex(i,j,k)-ex(i,j,k+1)-ez(i,j,k)+ez(i+1,j,k))*dt/(mu*jac(i,j,k))+chy(i,j,k)
             ! chz(i,j,k)=(ex(i,j+1,k)-ex(i,j,k)-ey(i+1,j,k)+ey(i,j,k))*dt/(mu*jac(i,j,k))+chz(i,j,k)
          enddo
       enddo
    enddo

    !write(*,*)'chz1= ',maxval(chz)

    ! to enforce div ch =0, from the bottom to top
    chz(:,:,nz2+1)=0


    continue
    do k=nz2,nz1,-1
       do j=ny1,ny2
          do i=nx1,nx2
             ! chz(i,j,k)=-chx(i,j,k)+chx(i+1,j,k)-chy(i,j,k)+chy(i,j+1,k)+chz(i,j,k+1)
             chz(i,j,k)=( &
                  -chx(i,j,k)*jac(i,j,k)+chx(i+1,j,k)*jac(i+1,j,k) &
                  -chy(i,j,k)*jac(i,j,k)+chy(i,j+1,k)*jac(i,j+1,k) &
                  +chz(i,j,k+1)*jac(i,j,k+1) &
                  )/jac(i,j,k)
          enddo
       enddo
    enddo

    !	write(*,*)'chz2= ',maxval(chz)


    ! call upward continuation, from chz(:,:,nz1) ==> chx(:,:,nz1-1),chy(:,:,nz1-1)

    !call upward()

  end subroutine migrate

end module mod_migrate
